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Abstract: 

The neuroimaging technique three-dimensional polarized light imaging (3D-PLI) provides a 
high-resolution reconstruction of nerve fibres in human post-mortem brains. The orientations 
of the fibres are derived from birefringence measurements of histological brain sections assum¬ 
ing that the nerve fibres - consisting of an axon and a surrounding myelin sheath - are uniaxial 
birefringent and that the measured optic axis is oriented in direction of the nerve fibres (macro¬ 
scopic model). Although experimental studies support this assumption, the molecular structure 
of the myelin sheath suggests that the birefringence of a nerve fibre can be described more pre¬ 
cisely by multiple optic axes oriented radially around the fibre axis (microscopic model). 

In this paper, we compare the use of the macroscopic and the microscopic model for simulating 
3D-PLI by means of the Jones matrix formalism. The simulations show that the macroscopic 
model ensures a reliable estimation of the fibre orientations as long as the polarimeter does 
not resolve structures smaller than the diameter of single fibres. In the case of fibre bundles, 
polarimeters with even higher resolutions can be used without losing reliability. When taking 
the myelin density into account, the derived fibre orientations are considerably improved. 

Keywords: polarized light imaging; nerve fibre architecture; optics; birefringence; 

Jones matrix calculus; computer simulation 


1. Introduction 


Unravelling the architecture and connectivity of nerve fibres in the human brain is one of the greatest challenges in 
neuroscience. Over the past years, several methods have been developed to reconstruct the human connectome (lj|3). 
The neuroimaging technique three-dimensional polarized light imaging (3D-PLI) has been employed to reconstruct 
the three-dimensional architecture of nerve fibres in human post-mortem brains with a resolution of a few micrometres 
@0- 3D-PLI enables the investigation of the pathways of long-range fibre bundles as well as single fibres and thus 
serves as a bridging technology between the macroscopic and the microscopic scale. 

The spatial orientations of the nerve fibres are derived by transmitting polarized light through histological brain 
sections in a polarimeter and measuring their birefringence. To relate the measured signal to the fibre orientation, an 
effective model of birefringence is used which assumes that the fibre density is constant over the whole brain section 0 
and that the measured optic axis indicates the predominant fibre orientation (5|6]]. This assumption is based on various 
experimental studies on white matter which show that the average birefringence of parallel nerve fibres is negatively 
uniaxial and that the measured optic axis is oriented along the length of the fibres [ 7||l0|. 

The majority of nerve fibres in the brain consist of an axon and a surrounding myelin sheath. The cytoplasm of 
the axon contains tubular polymers (microtubules) and neurofilaments running along the length of the axon |T]]|T2} . 
The myelin sheath is formed by oligodendrocytes (glial cells) which are spirally wrapped around the axon. The cell 
membranes are bimolecular layers consisting of lipid molecules and membrane proteins. The membrane proteins 
are embedded in the bilayer or attached to the membrane surface [ T3-T3J, whereas the lipid molecules are oriented 
radially to the fibre axis (T5]- 17]. The cell organelles of the axon and the protein framework of the myelin sheath lead 
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to a weak positive birefringence with respect to the longitudinal fibre axis [7||8|[T(),, 17|p0|. The anisotropic structure 
of the lipid molecules causes a positive birefringence with respect to the radial fibre axis |7| Emm- 

The effective model of uniaxial negative birefringence that is currently used in 3D-PLI seems reasonable for suf¬ 
ficiently low optical resolutions. However, it might no longer be valid if the anisotropic molecular structure of the 
nerve fibres is resolved. In this paper, we investigated the limitations of the effective model in terms of the optical 
resolution of the polarimeter using numerical simulations. The simulations were performed with a modified version 
of SimPLI [6], a simulation method that models the birefringence of the fibres with the Jones matrix calculus and 
allows data to be generated from synthetic fibre constellations that is comparable to experimental data. In order to 
study and understand the most dominant effects that generate the birefringence signals in 3D-PLI, the anisotropic 
molecular structure of the nerve fibres was described by a simplified birefringence model with radial optic axes (mi¬ 
croscopic model) and the effective model of uniaxial negative birefringence by a birefringence model with axial optic 
axes (macroscopic model). To investigate the limitations of the effective model, the transition between the microscopic 
and the macroscopic model was investigated depending on the optical resolution of the imaging system. 


2. Three-dimensional polarized light imaging (3D-PLI) 

The neuroimaging technique 3D-PLI determines the orientation of nerve fibres in post-mortem brains at the microme¬ 
tre scale. The principles of 3D-PLI have been explained in detail by Larsen et al. (22j and Axer et al. (4|[5). This 
section describes the measurement and data analysis procedures that are relevant for this study. 


2.1. Measurement 

To determine the orientation of the nerve fibres, a post-mortem brain - obtained from a body donor in accordance 
with ethical requirements - is fixed in buffered formaldehyde for several months, frozen and cut with a cryotome into 
histological sections of 70 pm, which are measured with a polarimeter. For the 3D-PLI measurement, two state-of-the- 
art polarimeters with different optical resolutions and sensitivities are employed: The large-area polarimeter (LAP) 
has a pixel size of 64 pm and is mainly used for single-shot images of whole human brain sections. The polarizing 
microscope (PM) has a pixel size of 1.33 pm (i. e. down to small axonal diameters), which enables complex fibre 
constellations to be disentangled. 

The LAP contains a pair of crossed linear polarizers and a quarter-wave retarder (with its fast axis adjusted at an 
angle of —45° with respect to the transmission axis of the first linear polarizer), see Fig. [T^. The employed light 
source emits incoherent, non-polarized, diffusive light with a peak wavelength of 525 nm. During the measurement, 
the polarizers and the quarter-wave retarder are rotated simultaneously around the stationary tissue sample. For each 
rotation angle p = 0°, 10°,..., 170°, the transmitted light intensity is recorded by a CCD camera so that a series of 18 
images is acquired. 

The imaging principle works as follows: The quarter-wave retarder transforms the linearly polarized light from the 
first polarizer into circularly polarized light. The birefringent brain tissue induces an additional phase shift so that 
the outgoing light is elliptically polarized. The fraction of light that then passes the second linear polarizer depends 
on the local orientation of the optic axis of the birefringent tissue, which is assumed to coincide with the local fibre 
orientation. 

The polarimetric set-up of the PM is slightly different to the set-up of the LAP (the order of the optical elements 
is reversed and only the first linear polarizer is rotatable). However, the imaging principle and the signal analysis are 
similar (4j so that the following considerations are only described for the LAP. 


2.2. Signal analysis 

The measured light intensity of an individual pixel describes a sinusoidal curve across the acquired image series, 
which depends on the orientation of the fibres within this pixel (see Fig. [lb). A physical description of the measured 
light intensity profile can be derived with the Jones matrix calculus |23||24] l, assuming that the light is coherent and 
completely polarized and that the optical elements are linear. For simplicity, the derivation is shown for a single pixel 
at a certain rotation angle p. 

In the Jones matrix calculus, all optical elements in the polarimeter are represented by Jones matrices (cf. Fig. GF 
The Jones matrices of the crossed linear polarizers are given by p5): 


Px = 
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The Jones matrix of a wave retarder that is rotated by an angle i/r in counterclockwise direction and induces along the 
fast axis a phase shift 8 between the two orthogonal components of the light wave is given by [ 25J: 
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In the experimental set-up, the fast axis of the quarter-wave retarder is rotated by —45° with respect to the axis of the 
first linear polarizer. Thus, the quarter-wave retarder can be described by the Jones matrix of a rotated wave retarder 
as given in Eq. (2.2) with a rotation angle of \jf = —45° and a phase shift of 8 = 90°: 


M V 4 = M 90 o(—45°) = 2= / 


(2.3) 


Under the assumption that the birefringence of the brain tissue can locally be described as negatively uniaxial with 
the optic axis indicating the predominant fibre direction (effective model), the brain tissue can locally be represented 
by a wave retarder that introduces a phase shift 8 along the fast axis (fibre axis). During the measurement, the two 
polarizers and the quarter-wave retarder are rotated simultaneously around the specimen stage in counterclockwise 
direction by a rotation angle p. For simplicity, the equivalent case is considered in which the brain tissue is rotated 
by an angle (—p) in counterclockwise direction while the other optical elements are fixed. Thus, the brain tissue can 
be described by the Jones matrix of a rotated wave retarder as given in Eq. (2 .2) with phase shift 8 and rotation angle 
iff = (p — p, where (p denotes the in-plane orientation of the optic axis: 
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When light with an electric field vector Eq passes through the 3D-PLI set-up, the resulting output beam with electric 
field vector Ej can be described by multiplication of the associated Jones matrices. As the Jones matrix calculus 
cannot be used to describe the non-polarized light emitted by the employed light source, the Jones vector E x — P x - Eo 
is used to describe the horizontally polarized light after the first linear polarizer (cf. Fig. |T^): 

Et = Py ‘ ^tissue * /4 ’ Ex- (2.5) 

Using It ~ \Et\ 2 , the transmitted light intensity is calculated, yielding a sinusoidal intensity profile (see Fig. [l])): 

MP) = ~y (l + sm( 2 (p-<p)) sin<5), ( 2 . 6 ) 

where It.q ~ \E*\ 2 corresponds to the transmitted light intensity averaged over all rotation angles (here referred to as 
transmittance) and | sin 8 | to the peak-to-peak amplitude of the normalized sinusoidal intensity profile (here referred 
to as retardation). The phase shift 5 is given by (see Appx. 0 : 

2k 9 

8 « —t An cos a , (2.7) 
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where A is the wavelength of the incident light, t the thickness of the brain section, An the local birefringence of the 
sample and a the local out-of-plane inclination angle of the fibre. Thus, the intensity profile in Eq. ( |2.6| ) is a direct 
measure of the spatial fibre orientation defined by the direction angle (p and the inclination angle a (see Fig. |TJc). 

In order to compute transmittance, direction and retardation, the intensity profile is fitted by means of a discrete 
harmonic Fourier analysis ji] 26]. The inclination angle a is calculated from the measured retardation | sin 8 | by 
rearranging Eq. ( |2.7| ). The direction and inclination angles are combined to a unit vector indicating the local fibre 
orientation in three dimensions. Putting all unit vectors of several adjacent brain sections together, a three-dimensional 
volume of vectors is created and the fibre tracts are reconstructed with streamline algorithms. 


3. Simulation of 3D-PLI using the Jones matrix formalism 
3.1. Simulation model 

3D-PEI derives the nerve fibre orientations based on the fact that the average birefringence of parallel fibres is neg¬ 
atively uniaxial (7}|l0| and assuming that the orientation of the measured optic axis corresponds to the local fibre 
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light source 


Figure 1. (a) Measurement set-up of 3D-PLI (for the LAP): The brain tissue is placed between a pair of crossed linear 

polarizers and a quarter-wave retarder, which are rotated simultaneously by 18 discrete rotation angles p. The transmitted 
light intensity is calculated with the Jones calculus, in which each optical element is represented by a Jones matrix (bold 
symbols), (b) The normalized transmitted light intensity h(p)/h,o describes a sinusoidal curve for each image pixel. The 
phase (p corresponds to the local fibre direction angle and the peak-to-peak amplitude |sin5| to the local fibre inclination angle, 
(c) The three-dimensional orientation of a fibre is defined by the direction angle <p and the inclination angle a. 


orientation. To investigate the limitations of this effective birefringence model, a straight single fibre and a hexago¬ 
nal bundle of straight parallel fibres were simulated and the birefringence of the fibres was modelled according to a 
microscopic and a macroscopic model for different optical resolutions of the simulated imaging system. 

Microscopic model: The microscopic model of birefringence considers the anisotropic molecular structure of a 
single nerve fibre. To investigate and understand the predominant effects generating the birefringence signals in 3D- 
PLI, a simplified model of birefringence was chosen for the simulations. As stated in Sec.[l] the average birefringence 
of parallel nerve fibres is negative with respect to the longitudinal fibre axis. Therefore, the positive birefringence of the 
axon and the myelin proteins is weak as compared to the birefringence effects of the myelin lipids |8|9|18|20|21| . Since 
the exact contribution of the different birefringence effects to the overall birefringence is unknown, the birefringence 
effects of the nerve fibres were modelled by considering only the anisotropic radial structure of the myelin sheaths: 
The fibres were simulated as hollow tubes (representing the myelin sheaths) with positive birefringence and radial 
optic axes (cf. lower Fig.[2j)). The axons were considered to be non-birefringent. 

Macroscopic model: To compare the simulation results of the microscopic model with the effective model of uniaxial 
negative birefringence, a macroscopic model of birefringence was defined. According to the assumptions made in the 
effective model, a single nerve fibre was simulated as negatively birefringent with axial optic axes oriented along the 
length of the fibre (cf. upper Fig. [2j)). Dohmen et al. © used this simulation model to investigate the effect of 
crossing fibre constellations. As this study concentrates on straight parallel fibres, the macroscopic model only serves 
as a reference for the effective model to verify the simulations of the microscopic model. To ensure a better comparison 
with the microscopic model, the fibres were simulated as hollow tubes (and not as solid cylinders as in ©)• 

3.2. Simulation method 

The basic idea of the simulation method is to model the birefringent myelin sheaths as series of linear optical retarder 
elements which are represented by Jones matrices. By defining the direction of the optic axes (radial/axial), both the 
microscopic and the macroscopic model can be simulated. The simulation approach is based on the simulation tool 
SimPLI developed by Dohmen et ah (6). For this study, the simulation tool was extended by the microscopic model 
and modified such that various fibre configurations with individual orientations, radii and myelin sheath thickness can 
be realized. 

The simulation tool is based on several assumptions and simplifications: First of all, the use of the Jones matrix 
calculus requires linear optical elements and perfect polarizers (i. e. the outgoing light is assumed to be completely 
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polarized). Another assumption is that the incident light can be described by parallel rays of light with straight optical 
pathways, i. e. the light is assumed to be non-diffusive and refraction, diffraction and scattering are neglected. For 
this study, a parallel and straight beam of light seems a reasonable approximation for the LAP because the imaging 
system has a small numerical aperture (the acceptance angle of the objective lens is less than 1°) so that the camera 
only captures light rays that are almost parallel to each other. 

The simulation consists of several steps: 

1.) Generation of synthetic nerve fibres in a three-dimensional volume: The nerve fibres are modelled as hollow 
tubes representing the myelin sheaths (see Fig. [2^). In order to approximate the geometry of the fibres, the simulation 
volume is discretized into small cubic volume elements (called voxels ), as indicated schematically by the grid in Fig. 

Et- 


2.) Generation of a three-dimensional vectorfield: For sufficiently small voxel sizes, the birefringence of the myelin 
sheaths can approximately be described by assigning each myelin voxel j a unit vector that indicates the direction of 
the optic axis ((pj, af) within the myelin sheath. In the macroscopic model, the vectors are oriented parallel to the fibre 
axis. In the microscopic model, the vectors are oriented radially to the fibre axis (see Fig.[2]3). 


3.) Generation of a synthetic 3D-PLI image series: In order to model the birefringence effect of the myelin sheaths, 
each myelin voxel is represented by the Jones matrix of a rotated wave retarder. The retarder axis is aligned with 
the optic axis within the myelin voxel (see Fig. [2]c). The synthetic 3D-PLI image series is calculated analogously to 
the derivation of the sinusoidal intensity profile as given in Eq. ( |2.5| ), with M t i ssue being replaced by the product of N 
matrices representing the myelin voxels along the optical path: 


Et =Py- (M N -M N ~i •••Ml)-M ^/4 


E x . 


(3.1) 


The matrix Mj = M§.((pj — p) is the Jones matrix of a rotated wave retarder as given in Eq. (2.2) and represents the 
j-th myelin voxel. The rotation angle depends on the in-plane direction angle (pj of the optic axis and the phase shift 
Sj on the out-of-plane inclination angle Oy. The Jones matrices of the linear polarizers and the quarter-wave retarder 
are given by Eqs. and p.3| ). For each rotation angle of the polarimeter (p = 0°, 10°,..., 170°), all Jones matrices 
along the optical path are multiplied (see Fig. [2]:), yielding a series of 18 synthetic 3D-PLI images with a sinusoidal 
intensity profile for each image pixel. 


3.3. Simulation parameters 

The choice of the simulation parameters was inspired by real experimental conditions. According to typical dimensions 
of large nerve fibres in human white matter | T6[[27| - [29) , the diameter and the myelin sheath thickness of the simulated 
nerve fibres were chosen to be 15 pm and 2.5 pm, respectively (see Fig. [3^). The fibres were generated in a simulation 
volume with dimensions v x y x z = 64 x 64 x 70 pm 3 , corresponding to the pixel size of the LAP and the thickness 
of the brain section. The simulation volume was discretized into cubic voxels with a side length of Ax s i m . In a 


preliminary study (see later, Sec. |4.1| ), the optimal voxel size was determined to be Ax s i m = 0.1 pm, which was used 
for all following simulations. Note that the dimensions are given in micrometres to meet the experimental conditions. 
As only relative length scales matter for the qualitative simulation results, the units could be chosen arbitrarily. 

Since measuring the birefringence of the micrometre-thick brain sections is impossible with the employed set¬ 
ups and literature values are not given for the currently used preparation technique, an upper limit for the bire¬ 
fringence of the myelin sheaths An was estimated: Under the assumption that a brain section that is completely 
filled with a homogeneous birefringent material with in-plane optic axis (a = 0°) induces a maximum possible re¬ 
tardation (| sin <51 = 1 <(=> 8 = 7t/ 2), the upper limit of the birefringence was calculated by rearranging Eq. (2.7): 
An = A/(4 1) = (525nm)/(4 • 70 pm) « 0.001875. Note that the choice of An only changes the overall magnitude 
of the retardation and does not affect the simulation results qualitatively. In the macroscopic (microscopic) model, the 
myelin voxels were simulated with axial (radial) optic axes and negative (positive) birefringence with respect to the 
optic axes. 

The wavelength of the incident light was chosen to correspond to the peak wavelength of the LAP (A = 525 nm). 
To study only the birefringence effect of the nerve fibres, the fibres were simulated without any absorption. 
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a) 


b) 





Figure 2. Simulation method: (a) Generation of synthetic nerve fibres in a three-dimensional volume, (b) Generation of a 
three-dimensional vector field according to the macroscopic model (axial optic axes) and the microscopic model (radial optic 
axes), (c) Generation of a synthetic 3D-PLI image series (illustrated for a large fibre in the microscopic model): The simulation 
volume is discretized into small volume elements (voxels). Each myelin voxel (grey) is represented by the Jones matrix of an 
optical retarder (Mj) whose axis is oriented in direction of the optic axes (arrows). The polarizing filters of the 3D-PLI set-up 
(see Fig. do are also represented by Jones matrices. For each rotation angle of the polarimeter, all Jones matrices along the 
optical path (highlighted column) are multiplied. 


3.4. Simulation of the optical resolution 

To investigate the effect of different optical resolutions on the measured 3D-PLI signal, the synthetic 3D-PLI image 
series were downsampled using the open-source image processing programme Fiji (30): To account for the limited 
optical resolution of the polarimeter, the image series were first convoluted with a two-dimensional Gaussian filter 
with a standard deviation <7. Then, the effect of the spatial discretisation of the CCD chip was modelled by resampling 
the resulting images with a sampling factor / s (average when downsizing). To determine realistic parameters for a 
and / s , the imaging properties of the LAP were considered as a point of reference (see Appx. 0 - 

Based on these considerations, the synthetic 3D-PLI image series were downsampled with different parameter sets 
(see Tab. [T}, yielding images with different pixel sizes Ax. The pixel size of the downsampled images was chosen 
such that a multiple of the pixel size corresponds to the side length of the simulation volume (Ax = 64 pm/n, with 
n = 4, 8, 16, 32). The standard deviation was calculated as a linear function of the pixel size (a = 0.714 Ax, see Appx. 
0 and the sampling factor was calculated by dividing the pixel size of the high-resolution image series by the pixel 
size of the downsampled image (/ s = Ax s i m /Ax = 0.1 pm /Ax). In the following, the optical resolution of the imaging 
system will be given in terms of the pixel size, which defines the set of downsampling parameters (cr and / s ) in Tab. 
[T] Note that the simulation results will not change qualitatively as long as the ratio between the fibre dimensions and 
the downsampling parameters remains the same. 


Ax [pm] 

a [|am] 

fs 

2.00 

1.43 

1/20 

4.00 

2.86 

1/40 

8.00 

5.71 

1/80 

16.00 

11.43 

1/160 


Table 1. Downsampling parameters (selected values): To obtain an image with pixel size Ax, a two-dimensional Gaussian 
filter with standard deviation a = 0.714 Ax and resampling with sampling factor / s = 0.1 pm/Ax are applied to the image. 
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3.5. Calculation of the retardation curve 


The determination of the inclination angle a is challenging for 3D-PLI because the peak-to-peak amplitude of the 
measured intensity profile (| sin 51) is highly sensitive to noise and - amongst others - influenced by the density of 
myelinated nerve fibres (see below). 

In the standard 3D-PLI analysis, the inclination angle is calculated from the measured intensity profile assuming 
that the brain tissue can locally be described by the effective model of uniaxial negative birefringence. In order to 
investigate whether the effective model can be used to extract the correct fibre inclinations, the retardation computed 
from Eq. \2.1) was compared to the retardation values derived from simulations using the macroscopic and the micro¬ 
scopic model (see Sec. HD- For that purpose, the retardation images were calculated for different fibre inclinations 
and different optical resolutions, respectively. For a better comparison between the retardation values of the single 
fibre and the fibre bundle, only the pixel in the centre of each (downsampled) retardation image was considered for 
evaluation. If pixels at other locations had been chosen, the retardation values of the single fibre would have been 
influenced by boundary effects that do not exist for the fibre bundle or real brain tissue which are completely filled 
with fibres. The retardation values from the centre of each downsampled retardation image were plotted against the 
corresponding inclination angle, yielding a retardation curve for each downsampling step. The retardation curves 
were compared to the normalized retardation curve of the effective model (cf. Eq. \2.1\ ), in the following referred to 
as theoretical curve : | sin 51 = |sin ((;r/2)cos 2 a) |. 

To be able to compare different retardation curves, the retardation was normalized for each curve with the maximum 
retardation value, respectively: 


sin 8 


n 8 

= sin I — - 


2 <5 ma 


(3.2) 


As only birefringent material (mainly myelin) is responsible for the phase shift in Eq. ( |2.7| ), t describes not the 
thickness of the whole brain section but rather the local myelin thickness t m , i. e. the combined thickness of myelin 
sheaths along the optical path. Due to the inhomogeneity of brain tissue, the local myelin density of a brain section is 
less than 100 %, i. e. the maximum possible retardation is | sin(5 a= o°, m ax) | < 1. If the inclination is calculated under 
the assumption that the maximum possible retardation equals 1, the inclination angle will be overestimated. In order 
to obtain a more precise estimation of the inclination angle, a so-called myelin density correction was applied to the 
downsampled retardation images: 

In the case of the macroscopic model, in which the optic axes within one fibre have the same orientations, 8 
scales linearly with t m . In the case of the microscopic model, in which the optic axes within one fibre have different 
inclination angles, the upper limit of 8 scales linearly with t m as long as the optic axes of neighbouring myelin voxels 
have similar orientations (see Appx. 0 Thus, the dependence on the myelin density can be eliminated to the greatest 
possible extent by multiplying the phase shift 8 with a correction factor (t/t m ): 


| sin (<5 corr ) | 



(3.3) 


In order to apply the myelin density correction to the downsampled retardation images, t m was replaced by the com¬ 
bined thickness of myelin voxels along the optical path (afte r ap plying the Gaussian filter and resampling). The 
resulting retardation images were normalized according to Eq. (3.2), yielding | sin(5 corr )|. 


4. Simulation results 

4.1. Comparison of analytical and numerical solution 

To estimate the accuracy of the simulation results for the microscopic model, a single fibre with radial optic axes 
and perpendicularly incident light (see Fig. [3])) was generated for different voxel sizes (Ax s i m ) and the numerically 
computed phase difference between extraordinary and ordinary wave (A<P num ) was compared to the analytical solution 

(A^ana)- 

Assuming that reflection and refraction effects can be neglected so that associated extraordinary and ordinary wave 
follow the same pathway, Bear and Schmidt derived an analytical expression for the phase difference (T8): 


. ^ 2tt 
A<P ana = 


471 A f 

— ro An I arccos 



— arccos 



(4.1) 


7 









a) 


b) 


i k 




Figure 3. (a) Dimensions of the simulated fibre (cross-sectional view), (b) Simulation model for comparison with the analyti¬ 
cal solution: A horizontal fibre is simulated with outer radius r\ = 1.5 pm, inner radius = 5 pm, and radial optic axes. The 
light is incident perpendicular to the fibre axis at distance tq. The electric field vector of the ordinary wave (E 0 ) is oriented 
parallel to the longitudinal axis of the fibre. The electric field vector of the extraordinary wave ( E e ) is oriented perpendicular 
to the fibre axis. 


where F is the optical path length difference between extraordinary and ordinary wave, r\ the radius of the whole 
nerve fibre (outer cylinder), r 2 the radius of the non-birefringent axon (inner cylinder) and ro the distance at which the 
light is incident perpendicular to the fibre axis (see Fig.j3j)). 

In order to compute A<F num , the propagation of ordinary and extraordinary wave were simulated separately: In 
the case of the ordinary wave, the light is polarized parallel to the longitudinal axis of the fibre. In the case of the 
extraordinary wave, the light is polarized perpendicular to the longitudinal axis of the fibre (see Fig. [3J)). The phase 
for both the ordi nary wave ( <P 0 ) and the extraordinary wave (<P e ) was calculated from the corresponding electric field 
vector Et in Eq. (3.1): 


<P - arctan 


( mjaiA 

\M\Et\)J 


(4.2) 


The numerically computed phase difference A<F num = <P e — <P 0 was evaluated at various distances 0 < ro < 5 pm away 
from the centre of the fibre and compared to the analytical solution given in Eq. ( |4.1| ), with r\ = 7.5 pm, = 5 pm, 
A = 525 nm and An = 0.001875 (cf. Sec. |3.3| ). In order to study the impact of the spatial discretisation on the accuracy 
of the numerical solution, the simulation was performed for various voxel sizes 1.50 pm > Ax s i m > 0.06 pm. As a 
measure of consistency between the numerical and the analytical solution, the relative phase difference was calculated: 
(A4> ana -A<J> num )/A4> ana . 

Figures^ and [4]) show the relative phase difference plotted against ro for various voxel sizes Ax s ; m . As can be seen, 
the numerical solution fluctuates around the analytical solution for voxel sizes of 0.5 pm and less. With smaller voxel 
sizes, the numerical solution approaches the analytical solution (indicated by the dashed black line). This behaviour 
is especially evident when considering the mean of the absolute relative phase difference for each voxel size (see Fig. 
U): For a voxel size of Ax s i m = 1.5 pm (corresponding to one tenth of the fibre diameter), the mean absolute relative 
phase difference is about 12 %. For Ax s i m = 0.5 pm, it is about 6 % and for Ax s i m = 0.06 pm, it is only 0.8 %. This 
demonstrates that the simulation tool produces correct results. 

As a good compromise between computation time and accuracy, all following fibre simulations were performed 
with a voxel size of Ax s i m = 0.1 pm (corresponding to 1/150 of the fibre diameter). For this voxel size, the relative 
phase difference is no more than 4 % (see Fig. |4]3) and the mean relative phase difference is about 1.3 % (see Fig. HP 


4.2. Simulation of a single fibre 

In a preliminary study, the limitations of the effective model of uniaxial negative birefringence were first studied for 
a straight single fibre. The fibre was simulated according to both the macroscopic and the microscopic model with 
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Figure 4. (a,b) Relative difference between analytically and numerically calculated phase difference (A0 a na and A0 num ) for 
various voxel sizes Ax s j m , evaluated at different distances ro away from the centre of the fibre. For reasons of clarity, the results 
are presented in two diagrams: (a) Ax s i m = 1.50-0.21 pm, (b) Ax s i m = 0.17-0.06 pm. The dashed black lines indicate the point 
at which the numerical values match the analytical solution, (c) Mean absolute relative phase difference plotted against the 
voxel size Ax s i m . The arrow indicates the voxel size (Ax s i m = 0.1 pm) that is chosen for the fibre simulations. 


different inclination angles (a = 0°, 10°, ..., 90°) and different optical resolutions. The dimensions of the fibre and 
the other simulation parameters were chosen as described in Sec.|3.3| The retardation curves were calculated from the 


downsampled retardation images (without/with myelin density correction) and normalized as described in Sec. [33 
An example of downsampled and corrected retardation images can be found in Appx. [D] 

Figure [5] shows the dimensions of the simulated single fibre and the corresponding retardation curves (continuous 
lines) for both simulation models and different optical resolutions (according to Tab. [I]). The theoretical retardation 
curve of the effective model is indicated by a dashed black line. In the case of the macroscopic model, the uncorrected 
retardation curves (see Fig. [5^) are already very similar to the theoretical curve for all investigated optical resolutions. 
After the myelin density correction (see Fig. [5]c), all retardation curves match the theoretical curve exactly, indepen¬ 
dently of the optical resolution. In the case of the microscopic model (see Figs. [5]) and[5]i), the retardation curves 
for a pixel size much smaller than the fibre diameter (Ax < 2 pm) are inverted as compared to the theoretical curve 
for a < 90°, i. e. the microscopic and the macroscopic model yield totally different results. For intermediate pixel 
sizes (2 pm < Ax < 8 pm), the retardation curves are non-monotonic, i. e. the assignment of the inclination angle is 
ambiguous. Finally, for pixel sizes larger than the fibre diameter (Ax =16 pm), the uncorrected retardation curve (see 
Fig. [5J3) is similar to the theoretical curve. After the myelin density correction (see Fig. [5ji ), the retardation curve 
matches the theoretical curve almost exactly. 


4.3. Simulation of a fibre bundle 

In brain tissue, nerve fibres are usually organised in hexagonal close-packed fibre bundles G3- In order to investigate 
the effect of fibre bundles on the 3D-PLI signal, a hexagonal bundle of straight parallel fibres with an inter-fibre 
spacing of 1 pm was simulated (see Fig. [6^). In order to obtain comparable results, the same dimensions and simulation 
parameters were chosen as for the single fibre. 

Figure [6] shows the normalized retardation curves for both simulation models and different optical resolutions (ac¬ 
cording to Tab.[l]). The (downsampled) retardation images that were used to compute the corrected retardation curves 
of the microscopic model are shown in Appx.|D| 
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Figure 5. (a-d) Normalized retardation curves of a straight single fibre simulated according to the macroscopic model (a,c) 
and the microscopic model (b,d) for different optical resolutions. Graphs (a,b) show the uncorrected retardation curves, graphs 
(c,d) show the retardation curves after the myelin density correction. For reasons of clarity, only selected graphs are shown. 
The legend indicates the pixel sizes of the retardation images from which the retardation curves have been calculated. The 
pixel size Av of the downsampled retardation images determines the parameters used for simulating the optical resolution (see 
Tab.[l|. For better comparison, Av is also given in terms of the fibre diameter (d = 15 pm), (e) Dimensions of the simulated 
single fibre. 
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Figure 6. (a-d) Normalized retardation curves of a hexagonal fibre bundle simulated according to the macroscopic model (a,c) 
and the microscopic model (b,d) for different optical resolutions. Graphs (a,b) show the uncorrected retardation curves, graphs 
(c,d) show the retardation curves after the myelin density correction. For reasons of clarity, only selected graphs are shown. 
The legend indicates the pixel sizes of the retardation images from which the retardation curves have been calculated. The 
pixel size Ax of the downsampled retardation images determines the parameters used for simulating the optical resolution (see 
Tab. [TJ. For better comparison, Ax is also given in terms of the fibre diameter (d = 15 pm), (e) Dimensions of the simulated 
fibre bundle. 
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In the case of the macroscopic model, the uncorrected retardation curves (see Fig. [6^) are very similar to the 
theoretical retardation curve of the effective model (dashed black line) for all investigated optical resolutions. As 
compared to the retardation curves of the single fibre (see Fig. [5^), the retardation curves of the fibre bundle are closer 
to the theoretical curve. After the myelin density correction (see Fig. [6]c), the curves are almost identical. In the 
case of the microscopic model, the uncorrected retardation curves (see Fig. [6J3) are also closer to the theoretical curve 
as compared to the uncorrected retardation curves of the single fibre (see Fig. [5j)). The myelin density correction 
(see Fig. [6ji) makes only a small difference, especially for low optical resolutions. For the simulated fibre bundle, 
the transition between the microscopic and the macroscopic model already occurs for pixel sizes larger than the fibre 
radius (Ax > 8 pm). 


Discussion 


In 3D-PLI, the fibre orientations are derived under the assumption that the brain tissue can (locally) be described 
as a homogeneous and uniaxial birefringent material with the optic axis indicating the predominant fibre direction. 
Furthermore, the density of myelinated fibres is assumed to be the same for the whole brain section. In this paper, 
the limitations of this effective birefringence model have been studied for the first time. For that purpose, a single 
fibre and a hexagonal fibre bundle (with diameters d) were simulated based on the Jones matrix calculus, employing a 
microscopic and a macroscopic model of birefringence and different optical resolutions (defined by the pixel size Ax 
as given in Tab.[I]). 

The transition between the two models is apparent when analysing the retardation curves: For high optical resolu¬ 
tions (Ax << d), the radial optic axes of the microscopic model are resolved. In this case, the optic axes are oriented 
perpendicular to the longitudinal fibre axis so that the retardation curves are inverted as compared to the macroscopic 
model and fibres with high inclination angles are interpreted as flat fibres. The zero retardation value for a = 90° is an 
artifact arising from the fact that the retardation is evaluated at the centre of the retardation image which - in the case 
of vertical fibres - contains no myelin (cf. Fig. [8j upper right corner). For intermediate optical resolutions (Ax < d), 
there is a transition zone between the microscopic and the macroscopic model so that an unambiguous assignment 
between retardation and inclination is not possible. For sufficiently low optical resolutions (single fibre: Ax > d; fibre 
bundle: Ax > d/2), the microscopic and the macroscopic model yield similar results (see Figs. [5ji and[6]l) so that the 
effective model of uniaxial negative birefringence can be used to compute the fibre inclinations. 

Thus, for the simulated fibre bundle (consisting of five fibre layers with d = 15 pm), the effective model can be 
used to interpret LAP measurements (Axlap = 64 pm > d/2), but not to interpret PM measurements (Axpm = 1 .33 pm 
< d/2). However, the diameters of the simulated fibres represent an upper estimate of typical fibre diameters in the 
human brain. The diameters of myelinated nerve fibres range from 0.3 to 15 pm |T6|27-29) and the majority of the 
fibres (e. g. 80 % in the corpus callosum [27 ]) have diameters of 1 pm or below so that the condition Axpm > d/2 is 
still fulfilled. In addition, fibre diameters much smaller than 15 pm implicate that the measured brain section (with 
thickness 70 pm) contains much more fibre layers than the simulated fibre bundle. A comparison between the simulated 
single fibre and the fibre bundle suggests that the more fibre layers are located along the optical path, the smaller is 
the minimum pixel size for which the effective model is still valid. To verify this hypothesis, the limitations of the 
effective model should also be studied in terms of the number of fibre layers along the optical path. However, a larger 
number of fibre layers also increases the probability that fibres with different spatial orientations are measured within 
the same volume, which poses a major challenge for 3D-PLI (6j. In future studies, the limitations of the effective 
model should therefore also be investigated for non-parallel fibre structures. 

The simulations have shown that - in regions with parallel fibre structures - the effective model of uniaxial negative 
birefringence is valid for the employed optical set-ups. For imaging systems with very high optical resolutions, the 
effective model needs to be reconsidered. Even if the optical resolution is too high to extract the correct fibre incli¬ 
nations, 3D-PLI remains a valuable neuroimaging technique as the image contrasts of transmittance and retardation 
still provide detailed structural information on the two-dimensional nerve fibre architecture in large histological brain 
sections. 

The effective model that is currently used for the data analysis in 3D-PLI does not only assume parallel fibre 
structures, but also a uniform myelin density. The simulations have shown that the retardation signal is considerably 
influenced by the myelin density, which impairs the reconstructed fibre orientations. It could be demonstrated that 
the estimation of the fibre inclination is considerably improved by the myelin density correction which incorporates 
the local myelin thickness of the examined tissue into the calculation of the inclination angle. While the correction 
has a large effect on the retardation curves of the single fibre, the effect is smaller for the fibre bundle which is 
much more homogeneous than the single fibre. Thus, the myelin density correction is especially useful for regions 
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with an inhomogeneous density of myelinated nerve fibres (e. g. for transition zones between white and grey matter). 
In the case of the microscopic model, the correction does not work as well as for the macroscopic model because 
the retardation also depends on the direction of the radially oriented optic axes in the myelin sheath, but it is still a 
considerable improvement. In order to incorporate the myelin density correction into the 3D-PLI signal analysis, the 
local myelin thickness t m of the sample needs to be determined. The intensity values of the transmittance image seem 
to be a good measure of the local myelin thickness in brain tissue [5J. 

The purpose of this study was to explore and understand the most dominant effects that generate the birefringence 
signals in 3D-PLI. To fully understand the physical processes behind 3D-PLI and to improve the interpretation of the 
reconstructed fibre orientations, a direct comparison between simulation and experiment is required. The long-term 
aim should be to develop a simulation tool of 3D-PLI that considers all relevant effects needed for reproducing the 
experimental results. To this end, the simulation model should be extended step by step and the relevant effects should 
be identified. 

Although the simulations show that the simplified microscopic model can already be used to explain the effective 
negative birefringence of parallel nerve fibres, future studies should include the positive birefringence of the axon and 
investigate how this modification changes the transition between the microscopic and the macroscopic model. 

So far, only straight and parallel fibres have been investigated. To provide more realistic fibre models, the fibres 
should be simulated with varying fibre diameters, myelin sheath thickness and spatial orientations. As fibres with 
different spatial orientations pose a major challenge for 3D-PLI [6J, future studies should focus on investigating 
inhomogeneous, non-parallel fibre structures. To enable a direct comparison with the experiment, the simulated fibre 
configurations should be based on experimentally determined fibre structures. 

In addition to a more realistic fibre model, the propagation of light should also be simulated more realistically. In 
this study, the incident light was described by a parallel beam of light. However, in the experiment, the employed 
light source emits diffusive light, i. e. the sample is illuminated by light with slightly different angles of incidence. 
As the measured birefringence signals depend on the angle between the light wave and the nerve fibres, a non-zero 
angle of incidence changes the retardation curves. For the LAP, which has a small numerical aperture, the effect can 
presumably be neglected. However, for systems with higher optical resolutions and higher numerical apertures, the 
effect of a Gaussian distribution of incident angles should be investigated further. 

Moreover, the simulations were based on the Jones matrix calculus which is only applicable to completely polarized 
and coherent light. As the light source emits incoherent light and the polarizers are not perfect, the Jones matrices 
should be replaced by Muller matrices GED which enable to study partially polarized and incoherent light. 

Finally, the assumption of a linear optical pathway is a great simplification. The refractive index of the myelin 
sheath is higher than the refractive indices of the inner axon and the surrounding tissue |J|[32][33j which will cause 
refraction/reflection at the interfaces and scattering of light. In future studies, the effects of refraction and scattering 
on the measured birefringence signal should be investigated in more detail. As the used simulation tool (SimPLI) is 
based on a matrix calculus, other simulation approaches will be required to investigate such non-linear pathways. 

6. Conclusion 

In this study, we laid a theoretical foundation for 3D-PLI. The effective model of uniaxial negative birefringence, 
which is currently used to compute the nerve fibre orientations from experimental data, has been validated for the first 
time. Using simulations based on the Jones matrix calculus, we have shown that the effective model can be used for the 
employed optical set-ups, i. e. as long as the polarimeter does not resolve structures smaller than the diameter of single 
nerve fibres. The developed Jones matrix formalism for simulating 3D-PLI has proven to be a powerful tool to gain a 
deeper theoretical understanding of the physical processes behind 3D-PLI and to better interpret the experimental data. 
The simulations enable not only to validate the computational model of the fibre reconstruction, but also to optimise 
the experimental set-up and the measurement method. 
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A. Derivation of the phase shift 


When polarized light passes through the birefringent brain section, it is split into an ordinary and an extraordinary wave 
which both experience different refractive indices. The refractive index n e that the extraordinary wave experiences 
when passing through the birefringent tissue under an angle 0 with respect to the optic axis, is given by (34): 
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n e (B) 2 


1 2 „ 1 • 2 „ 
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where n a is the ordinary refractive index and ///. = //,,(0 = 90°) the principal extraordinary refractive index of the brain 
tissue. 

The birefringence of biological tissue (An = he —n 0 = 10 -3 ...10 -2 [ 351) is small as compared to the values of the 
refractive indices n 0 and ///, (n = 1.3-1.5 |33|). Therefore, a Taylor expansion can be applied to the function 
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in An = 0: 
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The same expansion can be done for (l /n 2 0 — 1 /n e (0) 2 ) in (A n(Q) = n e (6) — 1). With these Taylor expansions, 

Eq. (A.2) can be written as: 

An(6) ~ An sin 2 0. (A.5) 

Choosing a coordinate system in which the light propagates in the z-direction and the brain tissue lies in the xy- 
plane, the optic axis (oriented in the direction of the nerve fibres) makes an angle 6 with the z-axis, i. e. the out-of-plane 
inclination angle of the fibre is a = 90° — 0. With this definition follows: An(0) ~ An cos 2 a. 

Thus, when the light passes through a brain section of thickness t, the extraordinary wave experiences a phase shift 
with respect to the ordinary wave which depends on the inclination angle of the optic axis: 

2 7i 2n 9 

S = —tAn(6) ~ -r-t Aftcos^ a. (A.6) 

A A 

This is the formula of the phase shift as given in Eq. (2.7) . 


B. Derivation of the downsampling parameters 


In previous measurements, the optical resolution of the LAP was investigated by employing a USAF test chart which 
contains line pairs (lp) with different spacings (36) . From the measured line intensity profiles, the Michelson contrast 
C was computed: 


C = 


7m ax An in 
Tnax H - Tnin 


(B.l) 


where 7 max corresponds to the mean intensity of the maxima and 7 m i n to the mean intensity of the (local) minima 
in the line intensity profile (cf. Fig. [7])). The largest number of line pairs per millimetre that can just be resolved 
(according to the Rayleigh criterion) was determined to be 5.66 lp/mm, which corresponds to a width per line pair of 
/lap = 176.7 pm and a contrast of Clap = 20.1 %. A width per line pair of 157.5 pm yields a Michelson contrast of 
11.2%. 
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According to these measurement results, a test image with three lines (pixel size: 0.1 pm) and a line width of 
4ap/2 ~ 88.4 pm was created, and the downsampling procedure (Gaussian filter and resampling) was applied to 
the test image (see Fig. |7^). The sampling factor was calculated by dividing the pixel size of the test image by the 
pixel size of the LAP: / s ,lap = 0.1 pm/64 pm. To reproduce the measured contrast of the line intensity profile (see 
Fig#), a Gaussian filter with a standard deviation of Opm = 45.7 pm was applied. To avoid boundary effects and 
ensure a symmetric line intensity profile, the dimensions of the image (1216 pm x 1216 pm) were chosen such that the 
downsampled image consists of an odd number of pixels (19 px x 19 px). 



Figure 7. (a) Downsampling of a test image (grey values: black = 0, white = 1): A Gaussian filter with standard deviation 
Olap and resampling with sampling factor / S; lap are applied to the test image, (b) Line profile of the downsampled test image: 
The determined contrast C = (/max — Anin) / (/max +/ m in) matches approximately the contrast Clap obtained from experimental 
measurements. 


Based on the determined parameters for the LAP (Axlap, Clap, / s ,lapX downsampling parameters for imaging 
systems with other optical resolutions (see Tab.[l]) were derived: The ratio between the pixel size of the LAP image and 
the determined standard deviation of the two-dimensional Gaussian filter is Olap/Axlap = 45.7 pm /64 pm « 0.714. 
Analogous measurements of the PM yield a similar ratio between pixel size and standard deviation J36| . Assuming 
that this ratio is the same for all simulated imaging systems, the standard deviation of the two-dimensional Gaussian 
filter was calculated from the pixel size Ax of the resulting downsampled image: 

a = 0.714Ax. (B.2) 

After applying the Gaussian filter, the synthetic image series (with pixel size Ax s i m ) was resampled with a sampling 
factor of 


= Ax sim = 0.1pm 
s Ax Ax 

yielding a downsampled image series with pixel size Ax. 


(B.3) 


C. Dependence of the phase shift on the local myelin thickness 

Each myelin voxel of a simulated nerve fibre is represented by the Jones matrix of a rotated wave retarder as defined in 
Eq. ( |2.2| ). Depending on what kind of model is used (macroscopic or microscopic), the retarder axis is either oriented 
parallel or radially to the fibre axis (see Fig.|2j)). 

In the macroscopic model, the optic axes of the myelin voxels are all oriented in the fibre direction (<p, a) so that 
the voxels can be described by the same Jones matrix j8) with phase shift 8 and (3 = (p — p. When the light 
propagates through N voxels of myelin, the multiplication of the N corresponding Jones matrices yields (using Eq. 
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(2.2) and R(P)R(~P) = I): 
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Thus, the N myelin voxels with thickness At (along the optical path) and phase shift 5 can be replaced by one myelin 
voxel with side length (N At M t m ) and phase shift: 

, J2jl 2 K n 271 o 

8 = N8 *=* ~z — Ah (N At ) cos 2 a = — Ant m cos 2 a. (C.2) 

A A 

In other words, the phase shift 5 (and for small 8 also the retardation |sin5|) scales linearly with the combined 
thickness of myelin voxels (N At), i. e. with the local myelin thickness t m . 

In the microscopic model, the optic axes of the myelin voxels along the optical path all have different orientations 
(cpj, ccj), see Fig. Hf If the optic axes of neighbouring myelin voxels have a similar direction ((p 2 — (pi <C 1 and 
0 l 2 — 0C\ <C 1), the multiplication of the N Jones matrices of the voxels can be simplified. For (p 2 — (p\ 1, one can 

define ft — ft = 7721 <C 1 and the multiplication of a pair of rotation matrices yields: 


p(_o n K ,f> \ _ f cos(P 2 ) sin(p 2 )\ fcos(Pi) -sin(fr)\ 

1 ^ ' K{Pl) ~\sm(p 2 ) cos(p 2 )J Vsin(j8i) cos(jSi) ) 

(cos(p 2 ~p l ) sin(jS 2 -j 8 i)\ / 1 rj 2 A 

\-sin(p 2 ~Pi) cos(P 2 —Pi)J ~ i?2i 1 / ’ 


using a first-order approximation in 7721 . 

The multiplication of two Jones matrices yields (with Mj = Mg. (ft)): 


M 2 -Mi 



Ip 


m) 


/ e i(5i+52)/2 ^^(52-50/2 

( r . J?21 g-i(52-5i)/2 e -i(5\+&i)/2 


e iSl / 2 0 \ 

0 g~ i 5 i/ 2 ) 

) *(-&)• 




The multiplication of four Jones matrices yields (ignoring terms in the order of rj? ): 




/ gH^l+%+^3+5*)/ 2 

V J ?" 


-i(5 1 +5 2 +5 3 +5 4 )/2 


where the elements of the secondary diagonal are given by: 

r j / = ^^1(54+ 53 +52-50/2 + ^2^(54+53-52-50/2 +I j 43 ei(a»-%-%-«i)/2 > 
rj" = -rj 21 g-i(5 4 +53+52-50/2 _ T?32e -i(5 4 +5 3 -5 2 -50/2 _ ^ 


(C.3) 


(C.4) 


(C.5) 


(C. 6 ) 

(C.7) 


If the number of myelin voxels (i. e. the number of matrices Mj) is small, the elements of the secondary diagonal in the 
resulting matrix can be neglected for rjjk <C 1. If the number of myelin voxels is large, the arguments of the exponential 
functions in the secondary diagonal will take all possible values and cancel each other for r\jk ~ r\i m V j.k.i.m. In 
both cases, the multiplication of N Jones matrices yields: 

M N -M N _ v ..M^R{p N ) ( Q e _ i{5i+ ... +5N)/2 ) R(-Pi). (C. 8 ) 
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Thus, the N myelin voxels with thickness At and phase shift Sj can be replaced by one myelin voxel with thickness 
(NAt = t m ) and phase shift: 


5' = ^ Sj -it-A n At ^ cos2 °Lj ^ (C.9) 

7=1 A 7=1 A 

given that the optic axes of neighbouring myelin voxels have similar directions. 

The analytical considerations have shown that the phase shifts of individual voxels add together in both simulation 
models: In the macroscopic model, the phase shift scales linearly with the local myelin thickness t m . In the microscopic 
model, this is only true for the upper limit of the phase shift. The dependence on the local myelin thickness is taken 
into account in the myelin density correction (see Sec. ig. 
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D. 


Retardation images of the fibre bundle 
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Figure 8. Retardation images of the hexagonal fibre bundle for selected fibre inclination angles a, simulated according to the 
microscopic model. The retardation values have been computed from (downsampled) image series with different pixel sizes 
Av (according to Tab.[lJ after the myelin density correction. The pixel sizes are indicated by a square in the left bottom corner 
of the retardation images. To enhance the image contrast, a different scale bar is used for each inclination angle (the minimum 
value of the high-resolution retardation image is encoded in black, the maximum value in white). 
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